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(54) Title: AUTOSEGMENTATION / AUTOCONTOURING SYSTEM AND METHOD 

(57) Abstract 

A system and method is disclosed for automatically 
computing contours representing the boundaries of ob- 
jects in three-dimensional tomographic images that may 
be formed by computed tomography ("CT"), magnetic res- 
onance imaging ("MRI"), positron emission tomography 
('PET"), single proton emission computed tomography 
("SPECT"), or other appropriate methods. The system 
and method begin with a sample region of the object's 
interior and the single region is expanded in a stepwise 
fashion. At each step, a contour maximally matching the 
region's current edge, local gray-level gradient maxima, 
and prior contour shapes is determined. Upon comple- 
tion of region expansion, the object contour is set to that 
step-contour having the maximum value of an objective 
function summing contributions from region edges, gra- 
dient edges, and prior shapes. Both the region expansion 
and the boundary contour determination are formulated 
such that there is a guaranteed average minimum error in 

the determination of the contours. This contour is represented as a parametric curve in which the contour size and shape are specified by 
the values of the parameters. These parameters are independant variables of the objective function. The parameters also are considered 
to be random variables capable of encoding a distribution of contour shapes, and by assuming a particular distribution, the contribution of 
shape constraints to the object function can be computed. The resulting contour corresponds to the set of parameters for which the objective 
function is a maximum. 
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Field Of The Invention 

10 The present invention relates to systems and methods for automatically 

inferring boundary contours of organs, tumors, prostheses, or other objects of medical 
interest from two-dimensional and three-dimensional images of the physical patient 
anatomy from computed tomography imaging, magnetic resonance imaging, or the 
like. 

15 

Background Of The Invention 

Clinical Imperatives and Issues 

Often in the course of clinical diagnosis, the patient's internal anatomy 
is imaged to determine the extent to which disease has progressed. The diseased tissue 
20 may be evidenced by some variance from normal anatomy or function. Several imag- 
ing modalities are commonly used to generate pictures (or images) of a patient's anat- 
omy and function suitable for diagnostic and radiotherapy treatment purposes, or for 
surgical planning. These include conventional X-ray plane film radiography; computed 
tomography (**CT") imaging, which also uses X-rays, magnetic resonance imaging 
25 ("MRI"), which produces images of internal anatomy and information about physio- 
1 logical function; and nuclear medicine imaging techniques, such as positron emission 
~ tomography ("PET") and single photon emission computed tomography ("SPECT"), 
which produce images with combined anatomic and physiologic or biochemical infor- 
mation. 

30 A common property shared by all the imaging modalities just men- 

tioned is that the images are digital. That is, the images are represented as regular 
arrays of numerical values, which represents a physical measurement produced by a 
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scanner. If these images are two-dimensional ("2-D"), the discrete picture elements are 
termed pixels. However, if the images are three-dimensional ("3-D"), the discrete vol- 
ume elements are termed voxels. For 3-D imaging modalities, single slices or sections 
are composed of pixels, but those same picture elements are equivalently termed vox- 
5 els when considering a set of stacked images as a volume of data. 

The digital images from 2-D or 3-D imaging modalities are substan- 
tially exact maps of the pictured anatomy, so that each pixel value represents a sample 
of a property at a location in the scanner's, and, therefore, patient's, coordinate system. 
Thus, the distances between pixel/voxel centers are proportional and have meaning in 
10 the sense of real physical spacing in the patient anatomy. Moreover, the relative posi- 
tioning of the numerical array represents the proportional spacing of objects in the 
patient anatomy. 

The numeric value of each pixel represents a sample of a property at 
that location. In CT images, for example, the numbers are a measure of relative X-ray 

15 absorbing power, so that spaces inside the lungs are usually pictured as dark (low CT 
number) while bone is generally bright (high CT number). 

CT imaging and MRI are two of the most frequently used imaging 
modalities because both provide detailed pictures of the internal anatomy of a patient. 
The instruments that employ these imaging techniques provide data that in appearance 

20 is 2-D or 3-D. However, the 3-D images, as stated, are a collection of 2-D samples, in 
this form of slices or sections, of the anatomy that have been combined to create a 3-D 
images. More specifically, to recreate the 3-D images from the 2-D image samples, the 
physician, scientist, or other skilled professional must recombine the 2-D image sam- 
ples (slices or sections) of the anatomic elements (organs, tumors, surgically-implanted 

25 prostheses, etc.) A common way to recombine 2-D image samples to form 3-D images 
is to manually draw individual contours on a contiguous set of 2-D image slices or sec- 
tions using computer graphics. Once these manually drawn contours are made, they 
are assembled to accurately construct 3-D representations of organs, tumors, and the 
like. The resulting 3-D reconstructions convey to the viewer the relative sizes, shapes, 

30 and mutual spatial relationships among the anatomic elements in the same anatomical 
scale as the original. 
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In the 2-D context of a slice or section, the individual anatomic ele- 
ments may be represented by contours coinciding with each object's boundaries. Alter- 
natively, in the 2-D context of a slice or section, anatomy elements may be represented 
by 2-D templates identical in size and shape to the object. 2-D templates are patterns of 
5 pixels all having the same value which represent a single region in an image. A repre- 
sentation by 2-D region-templates or by 2-D edge-contours are equivalent, since either 
representation can be readily computed from the other. 

As stated, 3-D reconstructions of patient anatomy are most often pre- 
pared using computer graphics by manually drawing the individual contours on a con- 

10 tiguous set 2-D image slices or sections and then combining them. This method is 
referred to as contouring. Contouring is very time-consuming and labor intensive. The 
time and labor necessary to use this method increases significantly with the number of 
image slices, and the number and sizes of the organs, tumors, etc. in the anatomical 
area of interest. The quality of the contouring and the later produced 3-D images, 

15 depend on the resolution and contrast of the 2-D images, and on the knowledge and 
judgment of the physician, scientist, or skilled professional performing the reconstruc- 
tion. 

Three-dimensional radiation therapy treatment planning ("RTTP") is a 
medical procedure that currently makes the greatest use of 3-D reconstructions. This is 
20 even despite the labor and time required to contour the organs and tumors to generate a 
useful plan. In fact, the largest fraction of the plan preparation time involves contour- 
ing. 

An example of a manually contoured CT image slice or section is 
shown in Figure i generally at 100. In Figure 1, the manually contoured organs are 
25 liver 102, spleen 104, left kidney 106, right kidney 108, and spinal cord 1 10. 

Figure 2, at 200, shows an example of a 3-D reconstruction that uses as 
an element the 2-D slice or section shown in Figure 1 at 100. The reconstruction in fig- 
ure 2 is composed of contours from a contiguous set of slices or sections. In Figure 2, 
the 3-D reconstruction of the liver is at 202, the spleen is at 204, the right kidney is at 
30 206, the left kidney is at 208, and the spinal cord is at 210. 

Another method that may be used for forming representations of 
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organs, tumors, and the like is the segmentation method. Segmentation is the identifi- 
cation of image objects as distinct regions or segments of an image. This method also 
may be used to generate 3-D reconstructions of a patient's anatomy. 

According to the segmentation method, a decision is made with regard 
5 to the image contents as to whether a given pixel of a 2~D slice or section belongs to a 
specific set of organs, tumors, lesions, or other objects known to exist in that slice or 
section. Therefore, given that both contouring and segmentation may equally be used 
to generate 3-D reconstructions, contouring and segmentation are taken to have the 
same meaning for description purposes herein. 

10 Pixel/ voxel values, which exist as intensities or gray levels, and their 

distributions across an image form a useful set of properties for segmentation. In a typ- 
ical slice or section, the edges of objects that are shown usually are associated with 
large value differences with nearby pixel values. Further, the interiors of discrete 
objects tend to have relatively constant values. As such, discrete objects exhibit distinct 

15 gray level textures in such a manner that adjoining objects or regions with different tex- 
tures appear to have visible boundaries between them. Each of these qualities of edge- 
ness and texture are associated with one or more computational methods that may be 
used to generate a numerical value for that property. As quantified, these properties can 
be used to make decisions about the segment identity of individual pixels. 

20 Prior Art Segmentation 

A number of autosegmentation methods have been proposed in the 
prior art. These prior art methods may be separated into two principal types: (1) semi- 
automated segmentation methods in which physicians, technicians, or skilled profes- 
sionals direct or provide some needed information which is used to produce detailed 

25 contours, and (2) fully automated segmentation methods in which a computer based 
program develops the segmentation without requiring any human intervention. These 
methods will be described in greater detail subsequently, 
a. Semi-Automated Segmentation Methods 

A known semi-automated method is reported in Cline, H.E., Lorensen, 

30 W.E., Kikinis, R., and Jolesz, E, Three-dimensional segmentation of MR! images of the 
head using probability and connectivity, J ournal of Computer Assisted Tomog raphy, 
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14:1037-1045, 1990. This method also is disclosed in whole, or in part, in U.S. Patent 
4,751,643 issued to Lorensen, et al. and U.S. Patent 4,791,567 issued to Cline, et al. 

According to this semi-automated segmentation method, the operator 
selects image pixels from distinct tissue in each of two or more multi-modal MRI 
5 image sets. The computer-based program then computes class conditional probability 
density parameters based on the assumption of a multivariate, Gaussian distribution 
data model. From seed pixels, the various tissues elements are grown outward by set- 
ting adjacent voxel labels according to a maximum likelihood decision rule. The label- 
ing process is carried out over the 3-D data of the multi-modal MRI image set, 
10 followed by smoothing of the segmented region boundaries to reduce noise. 

Another semi-automated approach is reported in Kohn, ML, Tanna, 
N.K., Herman, G.T., et al., Analysis of brain and cerebrospinal fluid volumes with MR 
imaging, Radiology. 178:115-122 (1991). According to this approach, manually- 
guided segmentation of the feature space leads to image-space segmentation used to 
15 measure brain and cerebrospinal fluid volumes in MRI images. 

Yet another semi-automated approach is reported in Hohne, K.H. and 
Hanson, W.A., Interactive 3D segmentation of MRI and CT volumes using morpholog- 
ical operations, Journal Qf Computer Agisted TpmQgr<ipfry, 16(2):285-294, 1992. Fol- 
lowing this approach, interactive gray level thresholding, and morphological erosion 
20 and dilation operations are used to more sharply define the apparent folds of the brain 
surface in 3-D graphics. 

A semi-automated segmentation method is reported in DeCarli, C, 
Maisog, J., Murphy, D.G.M., et al., Methodfor quantification of brain, ventricular, and 
subarachnoid CSF volumes from MRI images, Journal of Computer Assisted Tomogra- 
25 phy, 16(2):274-284 (1992). This approach is directed to segmentation of major brain 
regions, e.g., the cerebral cortical hemispheres, cerebellum, etc., based on an analysis 
of gray level histograms derived from manual samples of pixel values corresponding to 
respective brain tissue types. 

Another semi-automated segmentation method is disclosed in Neal, 
30 A.J., Sivewright, G., and Bentley, R., Technical note: Evaluation of a region growing 
algorithm for segmenting pelvic computed tomography images during radiotherapy 
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planning, British Journal of Radiolog y. 67:392-395 (1994). The disclosed method is a 
region-growing segmentation approach. According to this approach, pixels adjacent to 
pixels in an initial patch centered on a graphics cursor are added to the initial patch if 
the gray level of the pixel is within one standard deviation of the gray level mean of the 
5 initial patch. 

The semi-automated approaches that have been described all rely on 
two assumptions about objects in the images: 1) the pixel value means and variances 
are at least approximately constant throughout the object, which means that the object 
has statistical homogeneity, and 2) the neighboring objects have significantly different 

10 pixel value means, variances, or gray level distributions, which gives rise to visually- 
evident boundaries. However, problems can arise if these two assumptions are relied 
on too heavily because in both CT and MRI imaging in which anatomic structure con- 
touring is preformed, the two assumptions are frequently violated. Therefore, com- 
puter-based programs dependent on regional statistical homogeneity and/or high- 

15 contrast boundaries are likely to perform very poorly. 

b. Fully Automated Segmentation Methods 

Fully automated, computed segmentation has been reported only for 
limited anatomic locations and/or narrowly-defined imaging protocols. In fully auto- 
mated, computed segmentation system, an imaging modality is used to produce the 

20 original images. Modalities such as MRI are preferred because they produce images of 
soft tissue and display neighboring organs at higher contrast than X-ray based imaging. 
Further, MRI scanners can be set to produce images emphasizing proton density or dif- 
ferent relaxation phenomena. Further, multi-modal MRI, in principle, can provide 
more information for each voxel. 

25 The few fully automated, computed segmentation techniques that have 

been reported have been directed to the segmentation of the brain gray matter, white 
matter, and cerebrospinal fluid ("CSF") spaces using multi-modality MRI. These 
approaches use statistical pattern recognition methods to distinguish the various mate- 
rials. The main fully automated computed segmentation techniques using multi-modal 

30 MRI are disclosed and described in Bezdek, J.C., Hall, L.O., Clarke, LP, Review of MR 
image segmentation techniques using pattern recognition, Medical Phvsics . 20:1033- 
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1048 (1993); Fletcher, L.M., Barsotti, IB., and Hornak, J.P., A multispectral analysis 
of brain tissues, Magnetic Resonance in Medicine. 29:623-630 (1993); Vaidyanathan, 
M., Clarke, L.P., Velhuizen, R.R, et al., Comparison of supervised MRI segmentation 
methods for tumor volume determination during therapy, Magnetic Resonance Imag- 
5 ins, 13:719-728 (1995). 

A different strategy for fully automated, computed segmentation is to 
map a labeled atlas onto patient data by nonlinear transformations, referred to as warp- 
ing. This technique will produce local correspondences between the atlas and individ- 
ual patient anatomies despite inter-subject anatomic differences. An example of this 

10 strategy is reported in Miller MI, Christensen, G.E., Amit, Y., Grenander, U., Mathe- 
matical textbook ofdeformable neuroanatomies, Proceedings of the National Academy 
of Sciences USA . 90:1 1944-1 1948 (1993). The Miller et al. article describes a proce- 
dure in which an elastic model of the brain anatomy is driven by data-overlap probabil- 
ities to warp brain atlas images onto MRI slice or section images. In this case, 

15 segmentation occurs by associating the image voxels with atlas tissue-type labels. 

Another example of a fully automated, computed segmentation method 
is described in Staib, L.H. and Duncan, J.S., Boundary finding with parametrically 
deformable models, IEEE Transactions on Pattern Analysis and Machine Intelligence, 
14:1061-1075 (1992). According to this method, deformable models are used to seg- 

20 ment MRI image data. U.S. Patent 5,669,382 to Curwen et al., also describes this pro- 
cess. 

A further example of a fully automated, computed segmentation 
method is described in Chakraborty, A., Staib, L.H., Duncan, J.S., Deformable bound- 
ary finding in medical images by integrating gradient and region information, IEEE 

25 Transaction on Medical Imaging. 15:859-870 (1996). The method set forth in the 
Chakraborty et al. article is directed to acquiring images of the left ventricle at inter- 
vals in the cardiac cycle which are then used to construct 3-D models of the left ventri- 
cle to study cardiac function. 

The prior art discusses many other examples of deformable models that 

30 are used for segmentation. Many of these examples are discussed in the review of such 
technology in Mclnerney, T and Terzoploulos, D., Deformable models in medical 
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image analysis, in Mathemati cal Methods Biomedical Image Analysis, pp. 171-180 
(1996). 

Wells, W.M.III, Grimson, W.E.L., Kikinis, R., et al., Adoptive segmen- 
tation ofMRIdata, IEEE Transactions Medical Imaging. 15:429-442 (1996) discusses 
5 another fully automated, computed segmentation approach. According to Wells et al., 
there is simultaneous segmentation of the brain and CSF in MRI, and simultaneous 
correction of spatially-varying MR signal inhomogeneities using the knowledge of tis- 
sue intensity properties and intensity inhomogeneities. An expectation-maximization 
(EM) algorithm is used to compute the segmentation on the corrected intensities lead- 

10 ing to the classification of gray and white matter in quite distorted images. 

Some of the fully automated, computed segmentation methods that 
have just been described have experienced some success in very specific imaging 
modalities or anatomic settings, which have generally been the brain. The extent that 
these methods have found success over wider anatomic settings has depended on 

15 access to, and use of, supercomputers to compute solutions. Frequently, these methods 
must incorporate explicit anatomic knowledge, or knowledge of the image formation 
process. The ability of these methods to segment a number of different types of tissue 
in images having variable qualities, on single processor computers, has yet to be dem- 
onstrated. 

20 

Summary Of The Invention 

The present invention is an autocontouring/autosegmentation system 
and method for automatically computing contours representative of the boundaries of 
anatomical objects in two-dimensional ("2-D") and three-dimensional ("3-D") tomo- 

25 graphic images generated using computed tomography ("CP*) magnetic resonance 
imaging ("MRI"), positron emission tomography ("PET"), single photon emission 
computed tomography ("SPECT"), or other appropriate method. The system and 
method are useful for determining the boundary contours of an object in multiple sec- 
tions of a 3-D tomographic image in a novel way. The contours generated according to 

30 the system and method of the present invention optimally match (1) local image gray 
level gradients, (2) the edges of the segmented object, and (3) the prior contour shape 
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of the object at issue. These contours substantially conform to anatomic shapes with 
greater accuracy than hand-drawn contours. 

The system and method of the present invention provide an objective 
method for contour construction based on quantitative measurements of properties of 
5 the images themselves. More specifically, the system and method of the present inven- 
tion combines region, edge, and shape information to provide more accurate contours 
than methods using either region-growing or active contour approaches alone, such as 
reported in Chakraborty, A., Staib, L.H., Duncan, J.S., Deformable boundary finding in 
medical images by integrating gradient and region information, IEEE Transactions on 

10 Medical Imaging , 15:859-870 (1996). Further, the system and method of the present 
invention produce accurate contours in many anatomic sites, in images of variable 
quality, and require user interaction only to generate the first contour of an object. 

According to the present invention, the contours that are generated are 
closed parametric curves such that the (x t y) coordinates of points on the curve are 

15 themselves continuous functions of arc-length t and the set of real-valued parameters, 
p, as provided for in Expression (1): 

Uj) = (x(p,t),y(p 9 t)) (1) 



20 



25 



30 



where, 

(x,y) = the Cartesian coordinates of a point on a curve. 

t = the arc-length distance of a point (x,y) from an origin point. 

p - a set or vector of real- value parameters. 

x(p t t) = the x coordinate based on parameters p and arc-length distance r. 
y(p,t) = the v coordinate based on parameters p and arc-length t. 

The functions that define (x t y) will generate curves of arbitrary shape 
and size depending on the values of the parameters p. The parameters p serve as inde- 
pendent variables of an objective function which takes on greater or lesser values 
depending on the correlation of the computed, parametric contour and the actual object 
boundary in the image. 

The objective function that has been referred to is determined by the 
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sum of functions of the a posteriori conditional probabilities of the computed, para- 
metric contour given: (1) the quality of the match of the computed boundary with the 
perimeter of the interior region of the actual object, (2) the coincidence of the bound- 
ary with the local gray level gradient maxima, and (3) the similarity of the shapes of 
5 the estimated boundary with previously-determined, section boundaries of the actual 
object. 

The objective function is based on a Bayesian formulation to insure that 
the maximum a posteriori (MAP) result will be associated with the minimum average 
error in the contour determinations. The formulation of the objective function as a 

10 function of a posteriori probabilities follows what is set forth in Chakraborty, A., 

Staib, L.H., Duncan, J.S., Deformable boundary finding in medical images by integrat- 
ing gradient and region information, IEEE Transactions QnMMcal Imaging, 15:859- 
870 (1996) which thereby becomes part of the segmentation and method of the present 
invention. This maximized objective function is in the form M(pJ g J r ), where, p is the 

15 set of real-value parameters, I g is the gray-level gradient image, and I r is the region- 
classified image. 

In order to generate a contour according to the method of the present 
invention, the method initially requires a sample of the interior of the object to be con- 
toured, called the region-of-interest ("ROF'), which is obtained either by automated 

20 comparison with a neighboring section previously contoured by this method, or by 
manual input using interactive computer graphics. The system and method of the 
present invention then expands the ROI in a stepwise fashion by: (1) adding one layer 
of pixels where possible to the region's perimeter, (2) renumerating the perimeter, and 
(3) determining the set of parameters which maximally satisfy the three criteria above. 

^ Expansion continues in this manner to an appropriate stopping point. Once this is 

done, the object contour is defined by the set of parameters corresponding to the maxi- 
mum value of the objective function over all the expansion steps. 

The expansion or growing of the ROI according to the present invention 
is accomplished by testing image pixels outside the region, but adjacent to the region's 

30 perimeter. The decision whether a pixel being tested belongs inside or outside of the 
region is made by a supervised classifier decision method. This method computes the 
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values of discriminant functions (one each for each possible outcome) and chooses the 
outcome corresponding the highest valued discriminant function. The discriminant 
function is characterized by the assumption that the pixel gray level properties form 
multivariate Gaussian probability distributions. This supervised classifier decision 
5 method is developed de novo for each ROI for each section of each object. The values 
of the discriminant functions are set by the values of pixel gray-level properties and an 
assumed statistical model for those properties. 

The present invention will be disclosed in greater detail in the remain- 
der of the specification referring to the drawings. 

10 

Brief Description Of The Drawings 

Figure 1 shows a 2-D CT slice or section image through the abdomen, 
with certain organ contours labeled. 

Figure 2 shows a graphical representation of a 3-D reconstruction of the 
15 abdominal anatomy (which incorporates the 2-D CT slice or section image in Figure 
1), with certain organ contours labeled. 

Figure 3 is a graphical representation of a target pixel and the 3x3 set of 
pixels that influence the properties of the target pixel. 

Figures 4A, 4B, 4C, 4D, 4E, and 4F show a series of synthetic pictures 
20 in which the boundary (black) to be contoured is defined a priori by a fixed set of 
parameters p, and the computed contour (white) corresponds to the best, or maximum 
M(pJ gt I r ) objection function, for varying numbers of Fourier harmonics, N. 

Figures 5 A, 5B, 5C, and 5D show a series of synthetic pictures in which 
the boundary (black) to be contoured is defined a priori by a fixed set of parameters p> 
25 (as used in Figure 4), and to which increasing levels of Gaussian noise has been added. 

Figures 6 A, 6B, 6C, 6D, 6E, and 6F show a series of pictures of a sec- 
tion of the left kidney in which at Figure 6A an interior sample is inscribed (black 
polygon) by interactive graphics, and a contour is computed to fit the kidney boundary 
(dashed), and Figures 6B-6F show successive sections with significant variations in the 
30 shape details and interior texture of the left kidney. 
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A Detailed Description of the Invention 

The present invention is a system and method for automatically gener- 
ating accurate 2-D and 3-D contours using maximum a posteriori (MAP) autocontour- 
ing/autosegmentatton. According to the present invention, the boundary contours of a 
5 3-D tomographic reconstruction of an object may be determined from a contiguous set 
of 2-D slices or sections. These boundary contours optimally match the local image 
gray level gradients, the edges of the segmented object with accuracy, and the shape of 
the prior contour of the object. The system and method of the present invention effec- 
tively use region, edge, and prior shape information to provide accurate boundary con- 
I® tours in a novel way. 

The present invention generates the maximized objective function 
based on MAP regional growing, MAP boundary detection, and knowledge of the 
parametric representation of boundary shape. Each of these areas will now be dis- 
cussed. 

15 

Maximum a posteriori (MAP) Regional Growing 

Minimum Bayes risk is reported in Van Trees, H.L., Detection. Estima- 
te and Modul^tJmTheory, Pyt I, Wiley, New York, 1968. According to this refer- 
ence, minimum Bayes risk is a decision criterion which guarantees, on average, that 

20 the results of classification decisions will have the lowest expected loss. That is, if each 
possible decision outcome is assigned a loss (or cost), Minimum Bays risk decisions 
will have minimum average loss. 

A slightly simpler criterion for classification decision making is mini- 
mum Bayes error. Minimum Bayes error is reported in Duda, R.O. and Hart, P.E., Pat- 

25 tern Classification and Scene Analysis, Wiley, New York (1973); and Fukunaga, K., 
Introduction to Statistical Pattern Recognition, 2nd Ed., Academic Press . New York 
(1990). This criterion, which is based on cost = 0 for a correct classification decision 
and cost = / for an incorrect classification decision, guarantees that the average error 
of classification will be the lowest possible. As used herein, minimum Bayes error is 

30 the criterion for reclassifying an ROI-adjacent, outside pixel as an inside pixel. In 
addressing minimum Bayes error in greater detail, it will be directed to consideration 
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of all pixels in the inside and outside pixel classes and the computation of a set of gray- 
level properties or features. A statistically complete description of each feature, in each 
class, is given by the feature's probability density or histogram. Bayes rules provide 
the method to compare the features' probability descriptions across classes to make a 
5 minimum error decision. 

In order to better understand these decision rules for use in the present 
invention, the following is provided. In this example, X will represent a vector whose 
component values are the gray-level-derived properties, or features, for a single pixel 
whose class membership is to be tested. At least some of the features for this single 

10 pixel will normally have different values depending on its class, so the probability den- 
sities becomes the class conditional probability density p(X\i), which is to mean the 
probability of observing X given class /. In addition, the relative probability that class i 
is observed relative to other classes is the a priori probability P(i). Lastly, the probabil- 
ity that class i is the correct class for observation X is the a posteriori probability P(i\X) 

15 which is related to the other two probabilities mentioned by Bayes Rule as set forth in 
Expression (2): 



20 



25 



30 



P(i|X)= p(X\i)P(i)/L k p(X\k)P(k) 



where, 

P(i\X) = the a posteriori probability that class i in the correct class given 
vector X 

p(X\i) = the conditional probability of observing vector X given class i. 

P(i) = the a priori probability that class i is observed. 

p(X\k) = the conditional probability of observing vector X given class k. 

P(k) = the a priori probability that class k is observed. 

Z k p(X\k)P(k) = the sum of the products of the probability of observing vector 

X given class k and the a priori probability that class k is 
observed. 

P(i\X) = the a posteriori probability that class / is the correct class given 
vector X 
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After reviewing Expression (2) and understanding that it can be applied 
to two classes, the decision rule for two classes, such as classes i and j, will be to select 
the class associated with the largest a posteriori probability. This decision rule is 
embodied by Expression (3): 

if P{i\X) > P{j\X) , decide class i, else decide class y (3) 



where, 

10 P(i\X) = the a posteriori probability that class / is the correct class given 

vector X. 

P(j\X) = the a posteriori probability that class j is the correct class given 
vector X. 

However, if there are more than two classes, the choice will be the class corresponding 
15 to the maximum according to Expression (4): 

if / = arg max[P(y |X) ] , decide i (4) 

; 

where, 

20 

arg max[P(y \X)] = the value of the class label j corresponding to the largest 
j 

a posteriori probability P(jiX). 
The decision rule at Expression (4) represents the maximum a posted- 
25 ori probability rule. If prior probabilities are unknown, then the conditional probabili- 
ties, p(X\i), which are referred to as data likelihoods per Duda, R.O. and Hart, P.E., 
Pattern Classification and Scene Analysis, Wiley, New York (1973), may be combined 
in the likelihood ratio test of Expression (5): 

if P( X \ l \ > ^ decide class i, else decide class / (5) 
30 P(X\j) 
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where, 

p(XU) = the conditional probability of observing vector X given class L 
p(X\j) =s the conditional probability of observing vector X given class / 
i? = a threshold value greater than zero (0), which can accommodate 

5 changes in the estimates of prior probabilities, or give a mini- 

mum-error decision on known data (according to Van Trees, H.L., 
D^tMLjstim^tipn^nd Modulation Theory, Part I, Wiley, New 
York, 1968). 

If, however, the prior probabilities can be estimated, as in the present invention, the 
10 MAP probability rule for all classes, as set forth in Expression (3) is the appropriate 
rule to use. 

A review of Expression (2) indicates that it can be written in a more 
convenient form for use in the method and system of the present invention. From 
Expression (3), it is understood that p(i\X) p(X\i)P(i) because the denomina- 
15 tor in Expression (2) is a constant. As such, Expression (3) can be written as shown in 
Expression (6): 

if p(X\i)P(i) > p(X\j)P(j) , decide class i, else decide class j (6) 



20 



25 



30 



where, 

p(X\i) = The conditional probability of observing vector X given class /. 

P(i) = The a priori probability that class / is observed. 

p(X\j) = The conditional probability of observing vector X given class/ 

P(j) = The a priori probability that class j is observed. 

For reasons of computational convenience, discriminant functions, rep- 
resented by gj( j, of the product p(X\i)P(i) are often used. If the gtf ) are monotonically 
increasing, they may be substituted in Expression (6) without changing the functional 
nature of the decision rule. Thus, a pixel with vector X is assigned to class i according 
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to Expression (7): 

if g;(X) > gj(X), i* j , decide class /, else decide class j (7) 

where, 

gi(X) = A function of the a posteriori probability P(i\X) . 
gj(X) = A function of the a posteriori probability P(j\X). 
One such discriminant function is set forth in Expression (8): 
10 gj(X) = In p(X\i) + In P(i) (8) 

where, 

In = The natural logarithm of the conditional probability of vector X 

given class L 

1 S In F(i) = The natural logarithm of the a priori probability of class L 
gj(X) = A function of the a posteriori probability P(j\X) . 

Parametric forms for the class conditional probability densities p(X\i), 
p(X\j) can lead to expressions that are more convenient to evaluate. Therefore, with 
2Q regard to the present invention, there is an assumption that the pixel features are a mul- 
tivariate, Gaussian distribution that have a class conditional probability density accord- 
ing to Expression (8) (which follows Duda, R.O. and Hart, RE., Pattern Classification 
and Scene Analysis, Wiley, New York (1973)): 



25 



p(X\i) 



1 



d i 



exp[-(|)(X - Mfz. l (X - M,.)] (9) 



where, 

d = The number of feature -components in vector X. 

Mi = The vector of feature means for class L 

=s The matrix of covariances for the features in class L 

= The inverse covariance matrix for the features in class i. 
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E/| = The determinant of the covariance matrix. 

T = Transpose of the given vector. 

Now, if Expression (9) is substituted in Expression (8), and the natural 
logarithm of the density is taken, the result is the discriminant function at Expression 
5 (10): 

SiiX) = ^(X-Mfzj\x~M t )-^ln (2rc)-Q)ln + In (/>,.) 

(10) 



10 



30 



where, 

X - Vector X. 



Mi = The vector of feature means for class i. 

I<j = The matrix of covariances for the features in class i. 

= The inverse covariance matrix for the features in class /. 

15 IL^I = The determinant of the covariance matrix, 

T = Transpose of the given vector. 

J = The number of feature-components in vector X. 

P(i) = The a priori probability that class i is observed. 

An understanding of Expression (10) reveals that the ln(27t ) term can be omitted since 
2 ^ it is a constant for all classes. 

Keeping in mind Expression (10), the MAP test for two classes i and; 
according to at least one embodiment of the present invention, is the decision whether 
to classify a given pixel to the class of pixels belonging to the sample polygon, such as 
class i, or to the class of pixels outside the sample polygon, such as class j, based on 
the previously provided in Expression (7), which for convenience is provided here: 

if gj(X) > gj(X), i * j, decide class /, else decide class j (7) 



where, 

g ( (X) = A function of the a posteriori probability P(i\X). 
gj(X) = A function of the a posteriori probability P(j\X). 
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According to Expression (17) the feature vector X = {x h ...x n } has as its 
^-components the numeric values of several gray-level-derived measurements on the 
set of pixels in a neighborhood about each pixel. At least the use of the first order gray- 
level properties are reported in Pratt, W.K., Digital Image Processing, 2nd Ed., Wiley, 
5 New York (1991). The numeric values can include the mean, standard deviation, skew- 
ness, kurtosis, energy, entropy, and the range, but also other texture measures. How- 
ever, it is understood that there can be more gray-level-derived properties for which 
numeric values may be generated. 

The neighborhood that has been referred to preferably is a 3x3-pixel 
10 set, or the 8 pixels adjacent to the given pixel, minus any pixels not inside the ROI. 
From the values of the X-components, the mean feature vectors M, and the covariance 

matrices S f « can be computed and inserted into the Expression (10) for use in the deci- 
sion rule at Expression (7). 

A representative 3x3 pixel set is shown graphically in Figure 3, gener- 

15 

ally at 300. In Figure 3, the pixel of interest is center pixel 302. The pixels that influ- 
ence the numeric values of center pixel 302 are the 8 pixels adjacent to it. These are 
pixels 304, 306, 308, 310,312,314,316, and 318. These 8 pixels and the center pixel 
form the 3x3 pixel set except for pixels at, or adjacent to, the boundary of the polygon 
which fall outside the ROI. 

20 

Autocontouring/autosegmentation according to the present invention 
imposes three constraints on the probabilistic growth of the ROI to match to an image 
object boundary. The first is that the initial ROI polygon be entirely within the object 
being contoured. If the initial ROI polygon spans the boundaries between different 
anatomic structures, the present invention will not expand the boundary properly. The 

25 

second is that the ROI grows only be accreting former outside pixels. As such, no 
inside pixels are allowed to revert to outside at any point in the growth process. The 
third constraint is that any outside pixels that become completely surrounded by inside 
pixels are converted to be inside pixels. Thus, the resulting ROI perimeter will define a 
simply connected object. Given the foregoing, the application of the decision rules and 

30 

the constraints in the context of the present invention will now be described. 
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Maximum a posteriori (MAP) Boundary Detection: 
MAP Objective Function 

On average, the shapes of objects appearing in two or more serial-sec- 
tion images satisfy three properties: (1) the object boundaries coincide with local max- 
5 ima in the magnitude of the gray level gradient, (2) the boundary-enclosed regions 
have nearly homogenous textures, and (3) profiles of objects in a given section will be 
similar to those in the adjacent sections. Homogeneous textures are gray-level proper- 
ties that have constant means and variances across the whole region. The degree to 
which any trial contour coincides with local gradients and segmented region-edges, 
and agrees with shapes of prior contours, depends on the details of its shape. 

According to the present invention, computed contours are represented 
by continuous functions of shape parameters. Thus, the (x,y) coordinates of the points 
on the contour are themselves functions of the arc-length distance along the curve, t, 
along with the set (or vector) of parameters p f (x(p t t) y y(p,t)). If the total length of the 
contour is F, then 0 < / < T . The values of the parameters are given by equations for 
the Fourier elliptic representation which are described in detail below in the section 
entitled "Parametric Representation of Boundary Shape." 

The contour parameters also serve as independent variables for the 
objective function that (1) measures the region-edge and gradient-edge overlap, and (2) 
the similarity of the current contour with prior contours. Properly defined, this objec- 
tive function assumes a maximum value for the set of parameters p that correspond to 
the contour which most satisfies the three criteria stated at the beginning of this sec- 
tion. The objective function used in the system and method of the present invention is a 
function of the conditional probability reported in Chakraborty, A., Staib, L.H., Dun- 
can, J.S., Deformable boundary finding in medical images by integrating gradient and 
region information. IEEE Transactions on Medical Imaging . 15:859-870 (1996). This 
objective function is in the form of Expression (11): 

P{p\I r >l g ) < n > 

where, 
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p = a set or vector of real-value parameters. 

I r = the region-classified image, 

= the gray-level gradient image. 

The probability depicted by Expression (11) is the probability of 
5 obtaining the contour with parameter vector p given the region-classified image I n and 
the image of the scalar magnitude of the gray level gradient I g . The components of 
parameter vector p are treated as random variables, and parameter vector p is assumed 
to have a multivariate, Gaussian distribution probability density. 

It is reported in Chakraborty, A., Staib, L.H., and Duncan, J.S., Deform- 
^ able boundary finding in medical images by integrating gradient and region informa- 
tion, IEEE Transactions on Medical Imaging . 15-859-870 (1996) to use Bayes rule and 
the relationship between joint and conditional probabilities to derive the expression for 
P{p\I r I g ) in terms of measurable image quantities. The expression that is derived is 
15 at Expressions (12A, 12B, and 12C): 

- p ~w5t <12A) 



20 



25 



P { I r \I g ,p)P(p,I g ) (i2B) 



PU g , i r ) 



P(l r \l g ,p)p(p\j g )p(p) 



PU g ,I r ) (12C) 



where, 

P(p, I r , I g ) = The joint probability of contour vector p, and / ? . 

P{1 r \I g ,p) = The conditional probability of I r given I g and p. 

30 P(p, I g ) = The joint probability of contour vector p and f g . 

P{I g , l r ) = The joint probability of contour vector I g and I r 
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P(/?|/^ ) = The conditional probability of/? given 

- The a priori probability of the contour p. 
I r = The region-classified image produced by the MAP region growing 

j method described above. 

= The gray-level gradient image. 

In Expression (12A)~(12C), the denominator is constant, and further- 
more, the natural logarithm is a monotonically increasing function. Therefore, an 
equivalent expression for the boundary parameter values, which will be expressed as 
10 p* and maximize P(p\I r , I g ) , can be derived from Expression (13): 



p* = argmax[lnP(/ r |/ r p) + lnP{p\I g )+ InP(p)] 
p 



(13) 



15 



20 



where, 



argmax[ ] 
P 



PV r \I g >P) 
lnP(p\I s ) 



boundary parameter values. 

The argument corresponding to the largest term within the 
brackets. 

The conditional probability of I r given I g and p. 

The natural logarithm of the conditional probability of p given 



25 



30 



InP(p) = The natural logarithm of the a priori probability of values for 
P- 

Referring to Expression (13), the first term is the natural logarithm of 
the probability of a region image I r given gradient image I g and contour parameters p. 
The second term is the natural logarithm of the probability of obtaining a contour 
parameters p given the gradient image I g . The third term is the natural logarithm of the 
probability of a given contour parameters p. Obtaining actual probability measures is 
difficult because sufficient information about the space of possible outcomes is not 
available. Chakraborty, A., Staib, L.H., and Duncan, J.S., Deformuble boundary find- 



SUBSTITUTE SHEET (RULE 26) 



WO 00/08600 



-22- 



PCT/US99/11710 



ing in medical images by integrating gradient and region information, IEEE Transac- 
tions on Medical Imaging . 15:859-870 (1996) addressed this issue by rewriting 
Expression (14) to make explicit the dependencies of the terms on the various image 
properties, in a way related to, but not dependent on any probabilities. This resulted in 
5 Expression (14): 

arg maxM(/?, I g , I r ) = arg m&x[M pHor (p) + M gradien( (I g , p) + M region (I r p)] 
P P 

(14) 

10 where, 



15 



20 



25 



arg maxM(/?, I gy I r ) = the contour parameter vector p that corresponds to the 
P 

largest value of the objective function M of /?, I g , and I r 
= the value of the contour parameters p corresponding to 



arg max[ ] 
P 



M 



prior 



(P) 



M gradient^ 1 g » P) 



M region^ 1 \ • P) 



the largest value of the term in the brackets. 
= the function of the similarity of contour parameters p 

and the parameters of the corresponding contour of a 

neighboring section. 
= the function of the similarity of the gradient maxima and 

the contour specified by parameters p. 
= the function of the similarity of the classified region 

edge and the contour specified by parameters /?. 



According to Expression (14), the first (prior) term biases the boundary 
toward a particular distribution of shapes generated from prior experience. The second 
(gradient) term contributes the most when the parametric boundary /?, defined as the 
discrete boundary x(p,t),y(p,t) y most closely matches the coherent edge features in I g , 
30 The third (region) term is maximized when the parametric boundary p most closely 
matches the edges of the segmented region. The determination of these individual 
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terms will be discussed in detail below. 

a. Region Edge Contributions to the MAP Function: M region (p,I r ) 

The value of the M region (pj r ) terra in Expression (14) depends on the 
match of the parametric boundary with the edge of the region. Chakraborty, A., Staib, 
L.H., Duncan, J.S., Deformable boundary finding in medical images by integrating 
gradient and region information, IEEE Transactions on Medical Imaging. 15:859-870 
(1996) describes a method to maximize the exactness of the fit. This method rewards 
the boundary that contains as much of the inside region as possible, and penalizes the 
boundary that includes any of the outside pixels. The desired result is obtained by inte- 
grating over the area of a region template A p according to Expression (15): 



M„ gi0H {I„p) - JJ^ I{x,y)dA p (15) 

where, 

I r = the region-classified image. 

p s= a set of real-value parameters. 

A p = the area of the region template. 

In evaluating Expression (15), it is understood that the template pixels 
20 inside the region are set equal to +1 and outside pixels are -1. Thus, M region (pj r ) is 
maximized when p conforms to the edge of the region of +l's. Because the area inte- 
gral must be evaluated many times in the maximization of M(pJ g J r ), Chakraborty, A., 
Staib, L.H., Duncan, J.S., Deformable boundary finding in medical images by integrat- 
ing gradient and region information, IEEE Transactions on Medical Imaging . 15:859- 
25 870 (1996) describes an alternative integration method based on Green's Theorem 
which can be evaluated more rapidly. This alternative method will now be described. 

The alternative integration method based on Green's Theorem, results 
in the decomposition of the integral area into two line integrals which can be evaluated 
more rapidly. Green's Theorem, as reported in Kaplan, W., Advanced Calculus, Addi- 
30 $on- Weslev. Reading, MA, 1952, specifies that the area integral according to Expres- 



SUBSTITUTE SHEET (RULE 26) 



WO 00/08600 PCT/US99/1 1710 

-24- 



sion (16) can be written as the sum of two line integrals, as shown in Expression (16): 
lS A / r (x,y)dA= \\ c [* r (*o0§~ + M r (x#)&]dt (16) 

where, 

I r = the region-classified image. 

(x,y) = the x t y coordinates of any point in I r 
A p = the area of the region template. 

jq N r9 M r = the auxiliary functions defined in Expressions (17) and (18) 

below. 

t - the arc-length distance of a point (x,y) from an origin point. 

In view of Expression (16), the two terms to be integrated, M r and N r , 
may be written as Expressions (17) and (18): 

15 

M r (jc,y) = fl r {z,y)dz (17) 



N r (x,y) = -fl r (x,z)dz (18) 
20 J„ 

where, 

I r = The region-classified image. 

(x,y) = The x and v coordinates of a point. 
25 z = Dummy variable of integration. 

The image I r contains a template of the segmented region with pixel 
values of one (1) for pixels inside the region, and -1 for pixels outside the region. Eval- 
uating Expression (16) from the discrete image data is traversing the parametric 
boundary curve C p (in the plane olA p ) by summing values from 2-D arrays containing 

30 
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the sums M n N r as set forth in Expression (19): 

K 



10 



15 



1 ^ 

M regionVr*P) ~ 5 2 O*. )>(/>. O*) 



2 ^ ^r^i-a^vr-^ Af 
* = 1 

Ax(p y t) k 
M r (x(pj) k ,y(p, t) k ] 

(19) 

where, 

M region (h > P) - ^ e function of the similarity of the region-classified edge and 
the contour specified by parameters p. 

K = the length of the contour in pixels. 

N r , M r = the auxiliary functions defined in Expressions (17) and (18) 
above. 

x(p,t) = the x coordinate for specific values of parameters p and arc-length 
distance l 

y(p,t) = the y coordinate for specific values of parameters p and arc-length 
distance t. 

Ax = the finite difference computed from actual x values at pixel loca- 

tions. 

20 Ay = the finite difference computed from actual y values at pixel loca- 

tions. 

At = the finite difference computed from actual t values at pixel loca- 

tions. 

This sum set forth in Expression (19) is over all the /^-points in contour 
C p , and the A-differentials are evaluated by taking discrete differences. Discrete forms 



30 
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of M r , N r are given by Expressions (20) and (21), respectively: 



M r (x,y) = ]T / r (z»v) 

z = 0 



(20) 



Z = G 



(21) 



10 



where, 
N r ,M r 

fay) 



15 described 



= the auxiliary functions defined above. 

= the (x>y) coordinates of a point. 

= the region-classified image. 

- the index of sums over x or y coordinates. 

Noting the foregoing, the first term of Expression (14) has been 



Gradient Edge Contributions to the MAP Function 
The second term in Expression (14), M gradient (pJ g h depends on the 
coincidences of the parameterized boundary with edges in the image appearing as 
coherent features in the scalar gradient of the original image gray levels. The gradient 
20 term is a contour integral whose domain is the parametric contour C p realized by the 
discrete boundary [x(p,t)>y(p,t)]. In reviewing Staib, L.H. and Duncan, IS., Boundary 
finding with parametrically deformable models, TEEE Transactions on PattemAnalYi 
sis and Machine Intelligence * 14:1061-1075 (1992), it is found that it is assumed that/g 
contains a zero-mean noise process and independent boundary pixels making it possi- 
25 ble to evaluate M gradient (pJ g ) as the line integral according to Expression (22): 



M 



gradient 



(P,r g ) = J I g U{p,i),y{p,m 

G p 



dt 



(22) 



30 



where, 



= The gray-level gradient image. 
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10 



15 



20 



P 

x(p,t) 

y(p,t) 



a constant. 

a set of real- value parameters. 

the x coordinate for values of parameters p and arc-length distance 
t 

the y coordinate for values of parameters p and arc-length distance 
t. 



o = the noise variance. 

The likelihood of p representing the true boundary is proportional to the 
sum of the gradient values at all the points x(p,t),y(p,t) which will be used to evaluate 
the term M gradient (I gt p) over the K discrete contour pixels according the Expression 
(23): 



K 



M gradient^ r P) = £ 1 ' k^P > O*0>(/M)*) 
k = 1 



(23) 



where, 

'* 
K 

P 

x(p,t) 

y(p,t) 



25 described. 



= The gray-level gradient image. 
= the length of the contour in pixels. 
= a set of real- value parameters. 

= the x coordinate for values of parameters p and arc-length distance 
t. 

= the y coordinate for values of parameters p and arc-length dis- 
tance t. 

Given the foregoing, the second term of Expression (14) has been 



30 



Prior Contour Shape Contributions to the MAP Function 
The third term of Expression (14) is M prior (p). With regard to this term, 
it is assumed that the components of p = {p\>P2 P4n+i) torm a multivariate, Gauss- 
ian distribution. It is also assumed, however, that each component p t of the distribution 
set is statistically independent. The component probability densities Pip,) are set 
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according to Expression (24): 



P(Pi) = 



exp 



2a 2 i 



(24) 



where, 
Pi 



= the i-th component of the contour parametric vector p. 

=s the standard deviation for the component p i of the contour vector 



10 



m t - the mean value for the component p { of the contour vector p. 

Given Expression (24), the total probability for the multivariate, Gauss- 
ian distribution set is according to Expression (25): 



15 



^prior(P) = U P r(Pi) 
i = 1 



(25) 



The prior contour shape contribution to M(pJ g ,I r ) is found by combining Expressions 
(24) and (25) to get the products that is set forth in Expression (26): 

( (Pi-mf) 



20 



M priori) = Tl P riPi> JJ 



28 i J 



(26) 



25 



where, 
Pi 



- the i-th component of the contour vector p. 

= the standard deviation for the component p i of the contour vector 



= the mean value for the component p i of the contour vector p. 
Parametric Representation of Boundary Shape 

The functional form of the boundary parameterization is the Fourier 
elliptical representation. This is reported in Giardina, C.R. and Kuhl, F.P, Accuracy of 
curve approximation by harmonically related vectors with elliptical loci, Computer 
Graphics and Image Processing. 6:277-285 (1977); and Kuhl, FP and Giardina, C.R., 
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Elliptic Fourier features of a closed contour \ Computer Graphics and Image Process- 
ing , 18:236-258 (1982). According to these references, an object boundary contour is 
considered to be a closed, continuous curve V(p,t)> where t is the arc-length and p is the 
parameter vector. The curve has a total length T, such that 0 < t < T . The curve func- 
5 tion V depends on the discrete pixels located at (jc(p, t),y(p, t)) as set forth in the 
Expression (27): 



10 



15 



where, 

x(pj) = the x coordinate for values of parameter p and arc-length distance 
t. 

y(p,t) - the y coordinate for values of parameter p and arc-length distance 
t, 

p ~ a set of real- value parameters. 

t = the arc-length distance of a point (x(p,t),y(p,t)) from an origin 

point. 

The functions (x,(p,t)>y(p,t)) are periodic in arc-length distance, t> with 
2 q the period of the total curve length, Tt These are approximated by the finite Fourier 
series shown at Expressions (28) and (29): 

N 



x(p y t) = a 0 + X a « C0S "~T^ + ^" sm ^T^ 



n-1 



(28) 



25 " o 

0 = c* + 2, c^cos-^ + rf„sm— (29) 

In Expressions (29) and (30), the contour vector p is the set of Fourier 
coefficients {a 0 ,c 0 ,a i,b h c h d } ,...>a n> b w c n >d n }, and /V is the total number of Fourier har- 
monics. Maximization of the objective function, M(p,I g ,I r ), is carried out over the vec- 
tor {pi,P2> <P4N+2)* and the resulting contour is computed directly using Expressions 
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(28) and (29). The parameterization is global in that each parameter p t makes a contri- 
bution to (x(p,t),y(p,t)) at every value of t. 

As an example, to compute the contour parameters 
{a 0 ,c 0 ,a 1 ,b 1 ,cj,d I ,...,a n ,b n ,c n ,d„} from the coordinates (x(p,t),y(p.t)), Expressions (30)- 
^ (35) may be used 



15 



25 



30 



T * Ajc,r 2/jJt^ 2n7tr,._ ,-i 
2n 2 n 2 k = , At ^ T T J 



(30) 



10 T * Ajc,.r Inizt,. 2nnt k _, 

b n = — l —z y -A sin— ^ - sin 



„ 22^ a/ J r r J (3D 

2n tz k = i 



X a77 cos ^~ - cos-^- (32) 



2 2 ^ Af I r 



" = 2^V £ A*L " " "^""-1 
In it k - 1 K 

20 where, 

Ax = the finite difference computed from actual xffcj values at pixel 

locations. 

Ay = the finite difference computed from actual y(k) values at pixel 

locations. 

At = the finite difference computed from actual t values at pixel loca- 

tions. 

K = the length of the contour in pixels. 

T = the total contour length. 

t = the arc-length distance of a point x(k),y(k) from an origin point. 

The remaining terms, constant a 0 > c OJ are computed as follows accord- 
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ing to Expressions (34) and (35): 



20 



25 



30 



AT * * - 1 a .. * - 1 



(34) 



(35) 



where, 

10 Ax = the finite difference computed from actual x values at pixel locations. 

Ay = the finite difference computed from actual y values at pixel locations. 

At = the finite difference computed from actual t values at pixel locations. 

K = the length of the contour in pixels. 

T = the total contour length. 

15 t = the arc-length distance of a point x(k) f y(k) from an origin point. 

^ = 8, = 0. 

The set of Expressions (30) - (35) have been reported in Kuhi, F.P. and Giardina, C.R., 
Elliptic Fourier features of a closed contour, Computed Graphics and Image Processing, 
18:236-258 (1982). 

Given the Expressions (28) and (29), there is the need to know the number 
of Fourier harmonics, N, that are required for an acceptable approximation of a parametric 
contour for a given boundary. In this context, as the total length of the contour, T, 
increases, the number of Fourier harmonics possible increases with the support available 
for representation. However, there is an optimal number of Fourier harmonics for a given 
T. If the number of Fourier harmonics is increased from some small number, the corre- 
sponding contours display increasingly good agreement with the object boundary until the 
optional number is reached. However, there is a value of N beyond which the contours are 
increasingly degraded. 

The improved approximation occurs as N increases from a small number 
because the shape of the boundary is better captured by adding more information to the 
contour representation. Continuing to increase N past the optimum, however, inevitably 
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produces noisy contours because the Fourier coefficients are underdetermined with 
respect to the number of discrete pixels in the contour. This degradation is especially 
evident when the ratio of the number of contour-pixels K to the number of p-parame- 
ters (=4Ar+2) falls below 2.0. 
5 The behavior just described is consistent with the fact that small objects 

rarely have sharp corners needing high frequency representations by a large number of 
harmonics. Large objects such as the liver, for example, frequently have sharp turns in 
their boundaries but are well-represented by the large number of harmonics available 
because of the correspondingly larger supports (contour lengths). Thus, adapting the 

10 number of harmonics to the size of the object works because objects in CT and MRI 
are inherently (shape) band-limited, and, therefore, may be represented accurately by 
finite-length Fourier series. 

The relationship between the number of harmonics producing the best 
boundary-contour match and the contour length may be obtained, for example, by 

15 examining the contours produced for a series of abdominal organs with boundary 
lengths spanning an order of magnitude (see Figures 4A-4F). For each object, contours 
may be computed with increasing numbers of harmonics until the good matches are 
succeeded by noisy contours (see Figures 5A-5D). Plotting the number of harmonics 
needed to produce the best match versus the contour length results * in a log-linear 

20 relationship according to Expression (36): 

f 2 iflog 10 A:<21 

N*\ (36) 
\ int ( 12.0 * log 10 # - 22.4 ) if log io ii: > 2 1 



where, 

K = the length of the contour in pixels. 

int = the operation of taking the nearest integer of the argument in 

parentheses. 

N* = the number of Fourier harmonics. 

Expression (36) permits the system and method of the present invention to adapt to 
objects of varying size. 

The set of parameters p in Expression (1) serve as independent vari- 
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ables of the objective function M(p, I r , I g ), described above. The objective function 
assumes greater or lesser values (and significance) depending on the correlation of the 
computed, parametric contour and the actual object boundary of the image. The objec- 
tive function is based on the sum of functions of the a posteriori probabilities of the 
5 computed, parametric contour given: 

(1) the quality of the match of the computed boundary with the perimeter 
of the interior region of the actual object; 

(2) the coincidence of the computed boundary with local gray level gradi- 
ent maxima; and 

10 (3) the similarity of the shapes of the estimated boundary with previously- 

determined, section boundaries for the computed contour. 
The estimation and growing of the computed contour is based on the 
Baysian formulation of the objective function which insures that a maximum a posteri- 
ori (MAP) result will be associated with the minimum average error in computed, con- 
15 tour determinations. 

Now having described the components that result in the maximization 
of the computed contour, the implementation of the system and method of the present 
invention will be described. 

Implementation of the Preferred Embodiment 
20 The system and method of the present invention are used to generate 

computed contours that may be used for 2-D contouring and 3-D reconstructions, for 
example, of an anatomical areas of interest. The system and method of the present 
invention use data from the previously-contoured sections to generate contours for 
other sections, automatically, and without any user input. Generation of an initial con- 
25 tour does require interactive input, and that part of the system and method of the 

present invention will be described following the description of the automated contour- 
ing method. 

For the purpose of describing the implementation of the present inven- 
tion, some of the quantities that are generated during a current section's analysis take a 
30 new value during each step of the iterative region growing, and this will be indicated 
by the double superscript (k,l \ where, k is the section number and / is the iteration step. 
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A single superscript will indicate that the quantity is the final result of the contour 
determination for that section. 

The following is a description of an embodiment of the system and 
method of the present invention. 
5 The first action that takes place is that the present invention locates the 

pixels in the current section, I r (k * ! ~ 0) coincident with I r (k ' I} by classifying the current k~ 
section pixels with the (fc-1) texture classifier. The largest connected set of k pixels 
similar to the (kA) contour interior-pixels is the starting point for iterative region grow- 
ing. 

^ The present invention next makes a new texture classifier using the pix- 

els inside and outside of lj ktl ^ 0) . At that juncture, for each class /, the mean vectors 
m{ k) and inverse covariance matrices [Lf l Y k * are computed. Only one set of mean vec- 
tors and covariances are computed for each class i for the current section k. 

^ The present invention then computes the image of the gray-level gradi- 

ents for the current section I g (k K Once this is done, the present invention increments / 
by one, /=/+!. 

The present invention expands I r (ktl ~ ! * by one pixel layer to form I r (k,l \ 
and re-enumerates all the region perimeter pixels, 

20 

The present invention then forms from the l/ k,l) perimeter pixels, a set 

of contour parameters,//***^, and probability densities P(p (k * where * indicates that 
these form the initial value set before maximization of the objective function, and 
where the number of parameters (components of p) is varied with respect to the length 
25 of the contour to obtain the minimum-error match of contour to boundary. 

The present invention then maximizes the objective function in terms of 

the p parameters. The /-th cycle function is M* kfl) (p (k,l \ I g (k) J r (U) ) with the super- 
scripts indicating that parameters p (k>l * and region-increment I r (ktl > are specific to 
present iteration step. 

30 The present invention then tests the growth of the region to determine if 

it is at the maximum based on the decision rules for growth. This test is according to 
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the following procedures: 

If / > 2, and if / = the maximum allowed number of iterations or if 

< M* k,l ~ 2) y then the present invention will find the largest hS ktl) according to the 
Expression (37) using the values at the current interation step: 

5 

A# ( *V« / <» I r W ) « max[M ( ^(^ ;/) , / ( « /,<">)] (37) 

/ 

In Expression (37), the single superscript indicates that these are the 
optimal results to be saved and used for determining the (£+1) section. If, however, the 
predicate questions indicated immediately above are not satisfied, then the present 
invention will increment / by 1 and repeat the steps to this point in the growth of the 

region, beginning with the step in which I r (k ' I} is expanded by one pixel layer to form 
/ (U) 

15 Once the largest A4* kJ) is determined using Expression (37), the current 

it-th section results are saved. The information saved is 

(a) the template of interior of contour, I r (k) \ 

(b) the feature mean vector for each class i, m/ k} , /={ inside, outside}; 



10 



20 (c) the feature inverse covariance matrix for each class /, [Z i ] 

(d) the determinant of covariance matrix for each class /, 

(e) the shape parameters p (k) = {../?;.. } w and probability densities P(p <k) ). 
Now that this is completed, k is incremented by 1 for a new section and 

25 the process is repeated for this new section to be contoured. 

Having now described the method by which the ROI is grown to match 
the contour of a desired anatomical or other image, it is necessary to describe the 
method by which the initial ROI is generated. 

The initial contour is generated on the first section, k=0, by the user 
30 sampling the interior of the object to contoured. The ROI (or polygonal sample) that is 
generated is the starting point for the region growing that has been described which 
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uses a text classifier based on pixels inside and outside the initial ROL The texture clas- 
sifier has feature mean vectors and covariance matrices according to the following: 

(a) the feature mean vector for each class /, mf°\ /-{INSIDE, OUT- 
SIDE}; 

5 _} (o) 

(b) the feature inverse covariance matrix for each class /, ] 

(c) the determinant of the covariance matrix for each class i, 

[2 ( .] (0) . 

Each pixel outside of, and adjacent to, the ROI is tested, and if it is more 
like the inside pixels, its label is changed from OUTSIDE to INSIDE. The present inven- 
tion then re-enumerates all the perimeter pixels, and repeats the testing of all outside, 
ROI-adjacent pixels. This process is repeated until no further adjacent pixels pass the 
test. The last ROI is taken to be / r (0) . From the / r (0) perimeter pixels, a set of contour 

15 parameters and probability densities, P(/? ( * } ) are formed, where * indicates that 
these are initial values set before maximization of the objective function, and where the 
number of parameters (components of p) is varied with respect to the length of the con- 
tour to obtain the minimum-error match of the contour boundary. The image of the gray 
level gradient values is then computed, I g (Q). The set of contour parameters p(0) corre- 

20 sponding to the maximum of the objective function, Mf 0) (p ( *K I g (0 K I r (0) )> is according to 
Expression (38), 

MO) = argmaxrM< 0 V^V 0 > ( / r (0) )] (38) 
P ' 

25 

Now that the initial ROI is formed, the following parameters are saved so 
that the ROI may be grown using the region growing method starting with the next sec- 
tion, k=\: 

(a) the template of the interior of contour, //*" l) = / r (0) ; 
30 (b) the feature mean vector for each class /, m£ k ~ 1 } , mj (0) , /= { INSIDE, 

OUTSIDE}; 



SUBSTITUTE SHEET (RULE 26) 



WO 00/08600 



-37- 



PCT/US99/11710 



10 



(c) the feature inverse covariance matrix for each class i, 

(d) the determinant of the covariance matrix for each class /, 

P,f -" = Pf : 

(e) the shape parameters, p (k ' l) = { ../?,.. } (k " l) = { } (0) and the 

probability densities P(p (k ' l) ) = P(p m ). 
The system and method of the present invention have been described in 
detail above, now representative results on the operation of this system method will be 
described. 

Representative Results 

a, Computed contours as a function of the number of harmonics. 

Figures 4A-4F show a series of regions that have been grown using dif- 
^ ferent values of Fourier harmonics, N, More specifically, these figures consist of a 
series of synthetic pictures in which the boundary (black) to be contoured is defined a 
priori by a fixed set of parameters p and the computed contour (white) corresponds to 
the best or maximum M(pJ g J r ) objective function for the different numbers of Fourier 
harmonics, N. A comparison of the number of Fourier harmonics, N, and the RMS 
error values for Figures 4A-4F is shown in Table 1: 

Table 1 

Figure N RMS error 

4A 2 6.8 

4B 4 5.1 

4C 8 1.6 

4D 16 0.69 

4E 32 0.61 

4F 64 0.69 



20 



25 



30 



As shown in Table 1, the number of Fourier harmonics N~2 n y n=l to 6. 
Figure 4A shows boundary 402 and computed contour 404; Figure 4B shows boundary 
406 and computed contour 408; Figure 4C shows boundary 410 and computed contour 
412; Figure 4D shows boundary 414 and computed contour 416; Figure 4E shows 
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boundary 418 and computed contour 420; and Figure 4F shows boundary 422 and 
computed contour 424. A review of the match between the boundary and computed 
contour in the Figures 4A to 4F, respectively, reveals that there is improved matching 
as the number of Fourier harmonics is increased, with the best match (minimum error) 
5 occurring at 32 Fourier harmonics. 

In the set of matches defined by Figures 4A to 4F, it is readily seen that 
the match between boundary 402 and computed contour 404 in Figure 414 and bound- 
ary 406 and computed contour 408 in Figure 4B is not particularly good because much 
of the detail of the boundaries are not found in the contours. It is to be noted that Fig- 

10 ures 4A and 4B correspond to Fourier harmonic values of 2 and 4, respectively. It also 
is to be noted that the RMS (root mean square) error values are high for these harmon- 
ics. These high values are 6.8 for Figure 4A and 5.1 for Figure 4B. 

Figure 4C, which used 8 Fourier harmonics, has a significantly 
improved match between boundary 410 and computed contour 412. Table 1 also shows 

1 5 that the RMS error value has significantly decreased to 1 .6 compared to the RMS error 
values of 6.8 and 5 J for Fourier harmonics values of 2 and 4, respectively. Even with 
the improved matching in Figure 4C, some of the detail of the boundary is missing in 
the computed contour. 

Figures 4D, 4E, and 4F correspond to Fourier harmonics of 16, 32, and 

20 64, respectively. These Figures show significantly better matching between boundary 
414 and computed contour 416 in Figure 4D, boundary 418 and computed contour 420 
in Figure 4E, and boundary 422 and computed contour 424 in Figure 4F. In these suc- 
cessive figures, the computed contours that were grown with increasing numbers of 
Fourier harmonics become more defined compared to the boundaries. However, the 

25 incremental improvement is not overwhelmingly significant in evaluating the matches 
at Figures 4D, 4E and 4F. This also is indicated by the RMS values for Figures 4D to 
4F being about the same. As such, there is not significant improvement by using 16, 
32, or 64 Fourier harmonics although the match 4F, using 64 Fourier harmonics, visu- 
ally appears to be the best match. It is to be noted, however, that if the number of har- 

30 monies is increased above a certain level past the optimum, the result will in fact 
deteriorate in the match. 
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h. Computed contour accuracy versus noise 

Figures 5 A, 5B, 5C, and 5D show a series of regions that have been 
evaluated for noise. More specifically, Figure 5A to 5D show a series of synthetic pic- 
tures in which the boundary (black) to be contoured is defined a priori by a fixed set of 
5 parameters /?, which are identical to the set of parameters used in Figures 4A to 4F, and 
Gaussian distribution noise was added to the image. These images demonstrate the 
effect of noise on computed contour accuracy. Ten (10) Fourier harmonics are used for 
each of the figures. A comparison of the SNR (signal-to-noise ratio) and the RMS error 
for Figures 5A-5D is shown in Table 2: 
10 We 2 

Figure SNR RMS error 

5A <~ 1.2 

5B 2.0 L2 

5C 1.0 1.6 

5D 0.5 3.6 

15 

In Figure 5A, boundary 502 defines ROI 504. Initial polygon 506 is 
found in ROI 504 to be grown to boundary 502. The success in growing polygon 506 is 
influenced by the clarity of the boundary between the ROI 504 and the area 508 outside 
ROI 504. 

2Q According to Table 2, in Figure 5A the SNR is infinite. As is seen in 

reviewing Figure 5 A, there is a contrast between ROI 504 and area 508 outside bound- 
ary 502, but this high SNR also causes ROI 504 to have a significant gray-level value 
that is not desirable. This is true even though the RMS error is at a relatively low level. 

In Figure 5B, boundary 510 defines ROI 512. Initial polygon 514 is 
formed in ROI 512. In Figure 5B, the SNR is 2.0 and the RMS error has remained the 
same as it was when the SNR was infinite. In viewing Figure 5B, it is seen that ROI 
512 is in clear contrast to area 516 outside of ROI 512. This is a desirable situation to 
grow initial polygon 514 because this stark contrast will facilitate the growth to bound- 
ary 510 but not beyond it given the decision rules used by the present invention. 

In Figure 5C, boundary 520 defines ROI 502 and initial polygon 524 is 
formed in ROI 522. In Figure 5D, boundary 530 defines ROI 532 and initial polygon 
534 is formed in ROI 532. In Figure 5C there is not a great amount of contrast between 



25 



30 
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the gray-level values in ROI 522 and the area 526 outside of ROI 522. The same is true 
in Figure 5D, there is not a great amount of contrast between the gray-level values in 
ROI 532 and the area 536 outside of ROI 532. 

Figures 5C and 5D show situations in which the SNR is 1.0 and 0.5 
5 respectively. It also is seen that the RMS error values have increased slightly when the 
SNR has decreased to L0> but increased significantly when the SNR has decreased to 
0.5. These RMS error values are 1.6 to 3,6, respectively. In viewing Figures 5C and 
5D, it is evident that as the SNR is decreased from 2.0 to 1.0 and 0.5, the difference 
between the gray- level values of pixels inside the ROI and outside the ROI is much less 
10 distinguishable. Therefore, after considering Figure 5B, Figures 5C or 5D would not 
be desired environments for growing pplygon to match the ROI boundary. 

c« Automated contouring of the left kidney from CT images. 

Figures 6 A, 6B, 6C, 6D, 6E, and 6F show a series of pictures of a sec- 
tion of left kidney 602. In Figures 6 A to 6F, ten (10) Fourier harmonics were used. 
15 Figure 6A shows kidney 602 with interior inscribed (black polygon) 

604. Polygon 604 is inscribed using interactive graphics. From inscribed polygon 604, 
the present invention was used to compute the contour (dashed) 606 to fit the boundary 
kidney 602. 

Figure 6B to 6F show succeeding sections of left kidney 602 with sig- 
20 nificant variations in shape detail and interior textures. In each of these figures, the 
present invention was used to compute a contour that would fit the kidney section 
boundary. In Figure 6B contour 606 was computed, in Figure 6C contour 608 was 
computed, in Figure 6D contour 610 was computed; in Figure 6E contour 612 was 
computed; and in Figure 6F contour 614 was computed. 
25 The various kidney shapes and interior textures do not have any signifi- 

cant effect on the ability of the system and method of the present invention to compute 
the contours in Figures 6A-6F. 

The terms and expressions which are used herein are used as terms of 
expression and not of limitation. There is no intention in the use of such terms and 
30 expressions of excluding the equivalents of the features shown and described, or 

portions thereof, it being recognized that various modifications are possible in the scope 
of the present invention. 
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Claims: 

L A method for autocontouring a two-dimensional representation of a prede- 
termined bounded object that may be part of an image created by a predetermined 
modality, with the image being defined by a plurality of pixels having property values 
representative of what is on the image, comprising the steps of: 

(a) generating within the object an initial region of interest (ROI) with a bound- 
ary defined by a first supervised classifier based on properties of pixels inside the 
object compared with properties of pixels outside the object; 

(b) generating a second supervised classifier based on the properties of pixels 
inside the ROI; 

(c) expanding the ROI boundary using the second supervised classifier by eval- 
uating the properties of a layer of pixels adjacent to a current exterior boundary of the 
ROI based on comparing pixel properties of each pixel in the current exterior layer to 
the properties of the pixels within the ROI, reclassifying the pixels in the current exte- 
rior layer as pixels to be included as part of the pixels within the ROI if the properties 
of discrete pixels in the current exterior layer substantially match in a predetermined 
manner the properties of the pixels within the ROI, and re-numerating the boundary of 
the ROI based on any reclassified pixels; 

(d) generating a parametric representation of a contour based at least in part on 
two-dimensional coordinates of a predetermined number of points on the ROI bound- 
ary determined in step (c); 

(e) generating an optimal objective function and a corresponding parametric 
contour; 

(f) repeating steps (c)-(e) one exterior pixel layer at a time until a maximum 
match is captured; 

(g) selecting a largest-valued optimal objective function that has been gener- 
ated according to steps (e) and (0; and 

(h) creating a contour based of the objective function selected at step (g). 

2. The method as recited in claim 1, wherein the predetermined modality gen- 
erates contours in tomographic images. 

3. The method as recited in claim 2, wherein the predetermined modality 
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includes X-ray computed tomography (CT). 

4. The method as recited in claim 2, wherein the predetermined modality 
includes magnetic resonance imaging (MRI). 

5. The method as recited in claim 1, wherein the parametric representation is in 
a form of a finite Fourier series of two-dimensional coordinates that define the bound- 
ary of the ROL 

6. The method as recited in claim 5, wherein the finite Fourier series includes 
the use of a predetermined number of Fourier harmonics to minimize error between 
points the contour and points on the boundary of the ROI. 

7. The method as recited in claim 6, wherein the number of Fourier harmonics 
relates to a length of the boundary. 

8. A method for autocontouring a two-dimensional representation of a prede- 
termined bounded object that may be part of an image created by a predetermined 
modality, with the image being defined by a plurality of pixels having property values 
representative of what is on the images, comprising the steps of: 

(a) generating within the object an initial region of interest (ROI) with a bound- 
ary; 

(b) generating a supervised classifier based on predetermined properties of pix- 
els inside of the ROI and properties of pixels outside of the ROI for evaluating if pixels 
in a layer ad jacent to a current exterior boundary of the ROI should be included as part 
of an interior of the ROI; 

(c) expanding the boundary of the ROI using the supervised classifier by itera- 
tively evaluating exterior layers of pixels adjacent to the boundary of the ROI until 
there is maximum matching of the ROI boundary and the object boundary; and 

(d) generating a contour for the object according to the maximum expanded 
ROI boundary based on a maximization of an objective function. 

9. The method as recited in claim 8, wherein the predetermined modality gen- 
erates contours in tomographic images. 

10. The method as recited in claim 9, wherein predetermined modality includes 
X-ray computed tomography (CT). 

11. The method as recited in claim 9, wherein predetermined modality includes 
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magnetic resonance imaging (MRI). 

12. A method for autocontouring a three-dimensional representation of a prede- 
termined bounded object that may be part of an image created a series of two-dimen- 
sional images formed by a predetermined modality, with each of the two-dimensional 
images being defined by a plurality of pixels having property values representative of 
what is on the image, comprising the steps of: 

(a) generating within an object on a single image an initial region of interest 
(ROI) with a boundary defined by a first supervised classifier based on properties of 
pixels inside the object on the single image compared with properties of pixels of with 
the same object in an adjacent image; 

(b) generating a second supervised classifier based on the properties of pixels 
inside the ROI; 

(c) expanding the ROI boundary using the second supervised classifier by eval- 
uating the properties of a layer of pixels adjacent to a current exterior boundary of the 
ROI based on comparing pixel properties of each pixel of the current exterior layer to 
the properties of the pixels within the ROI, reclassifying the pixels in the current exte- 
rior layer as pixels to be included as part of the pixels within the ROI if the properties 
of discrete pixels in the current exterior layer substantially match in a predetermined 
manner the properties of the pixels within the ROI, and re-numerating the boundary of 
the ROI based on any reclassified pixels; 

(d) generating a parametric representation of a contour based at least in part on 
two-dimensional coordinates of a predetermined number of points on the ROI bound- 
ary determined in step (c); 

(e) generating an optimal objective function and a corresponding parametric 
contour; 

(f) repeating steps (c)-(e) one exterior pixel layer at a time until a maximum 
match is captured; 

(g) selecting a largest-valued optimal objective function that has been gener- 
ated according to steps (e) and (f); 

(h) saving the objective function selected at step (g); 

(i) repeating steps (a) - (h) for each two-dimensional image defining the three- 
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dimensional representation of the object; and 

(j) creating a three-dimensional contour based of the objective functions 
selected and saved at steps (g) and (h). 

13. The method as recited in claim 12, wherein the predetermined modality 
generates contours in tomographic images. 

14. The method as recited in claim 13, wherein the predetermined modality 
includes X-ray computed tomography (CT). 

15. The method as recited in claim 13, wherein the predetermined modality 
includes magnetic resonance imaging (MRI). 

16. The method as recited in claim 12, wherein the parametric representation is 
in a form of a finite Fourier series of two-dimensional coordinates that define the 
boundary of the ROL 

17. The method as recited in claim 16, wherein the finite Fourier series includes 
the use of a predetermined number of Fourier harmonics to minimize error between 
points the contour and points on the boundary of the ROL 

18. The method as recited in claim 17, wherein the number of Fourier harmon- 
ics relates to a length of the boundary. 

19. A method for autocontouring a three-dimensional representation of a prede- 
termined bounded object that may be part of an image created a series of two-dimen- 
sional images formed by a predetermined modality, with each of the two-dimensional 
images being defined by a plurality of pixels having property values representative of 
what is on the image, comprising the steps of: 

(a) generating within the object on a single image an initial region of interest 
(ROI) with a boundary; 

(b) generating a supervised classifier based on predetermined properties of pix- 
els inside of the ROI and properties of pixels outside of the ROI for evaluating if pixels 
in a layer adjacent to a current exterior boundary of the ROI should be included as part 
of an interior of the ROI; 

(c) expanding the boundary of the ROI using the supervised classifier by itera- 
tively evaluating exterior layers of pixels adjacent to the boundary of the ROI until 
there is maximum matching of the ROI boundary and the object boundary; 
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(d) generating an optimal objective function for the maximum matching ROI 
boundary; 

(e) saving the optimal objective function generated at step (d); 

(f) repeating steps (c) - (e) for each two-dimensional image defining the three- 
dimensional representation of the object; and 

(j) generating a three-dimensional set of contours based of the objective func- 
tions generated and saved at steps (d) and (e). 

20. The method as recited in claim 19, wherein the predetermined modality 
generates contours in tomographic images. 

21. The method as recited in claim 20, wherein the predetermined modality 
includes X-ray computed tomography (CT). 

22. The method as recited in claim 9, wherein the predetermined modality 
includes magnetic resonance imaging (MRI). 
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